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We  show,  via  scaling  arguments,  that  a  first  order  transition  occurs 
when  the  appropriate  susceptibility  (second  derivative  of  the  free  energy) 
diverges  proportional  to  the  volume  of  the  system.  This  result  holds 
rigorously  as  an  equality  or  an  upper  bound  for  a  very  large  class  of 
Hamiltonians.  The  divergence  is  maximal;  therefore  it  provides  a  useful 
way  of  locating  first  order  transitions.  Our  method  gives  accurate 
results  for  the  two  dimensional  Ising  model. 
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We  formulate,  by  means  of  a  scaling  argument,  a  general  criterion 
for  first  order  phase  transitions  In  terms  of  the  behavior  of  systems  of 
finite  size.  This  work  Is  a  natural  extension  of  the  Ideas  of  Nlenhuls 
and  Nauenbergl  on  first  order  transitions  In  Infinite  systems  and  those 
of  Nightingale^  on  calculating  the. critical  temperature  -( and  critical  - 
exponents)  for  higher  order  transitions  from  finite  size  results.  In 
common  with  the  latter  and  other  recent  works^  we  find  that  useful 
information  about  the  thermodynamic  limit  can  be  obtained  by  comparing 
finite  systems  of  different  sizes.  Further,  we  point  out  that  our  result 
is  rigorously  true  for  a  very  large  class  of  ferromagnetic  Ising  models^. 

Our  criterion  provides  a  way  of  estimating  the  temperature  at  which 
a  first  order  transition  occurs  from  finite  size  results.  It  also 
implies  that  at  a  first  order  transition  a  second  derivative  of  the  free 
energy  (such  as  the  magnetic  susceptibility)  diverges  proportional  to  the 
volume  of  the  system.  At  a  higher  order  transition,  the  corresponding 
argument  indicates  a  less  rapid  divergence. 

Our  criterion  is  exact  for  the  two  dimensional  Ising  model5.  We 
find,  from  finite  size  calculations,  that  the  convergence  Is  rapid  and 
one  can  accurately  locate  the  line  of  first-order  transitions. 

Our  argument  is  quite  general.  However,  for  purposes  of  illustra¬ 
tion,  consider  a  magnetic  sy-’tem  (e.g.  the  Ising  model)  of  N  spins  In 
dimension  d  at  temperatur*.  i  an  external  magnetic  field  h.  Let  (T,  h) 
be  the  total  free  energy,  f  =  F^/N  the  free  energy  per  spin,  and 
m(T,h)  the  magnetization  per  spin. 
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Now  consider  a  similar  magnetic  system  with  N1  spins  at  the  same  tempera¬ 
ture  T  but  some  different  external  field  h'.  Then 

FN  (T.h)  -  FN,  (T.h*)  +  6(T,h)  (1) 

If  h  »  h',  1g  Is  the  free  energy  difference  between  the  two  systems.  For 
h  f  h',  Eq  1(1 )  resembles  the  Inhomogeneous  scaling  equation  for  the  free 
energy  In  renormalization  group  theory®.  Dividing  by  N  we  have 

f(T,h)»4f(T,h*)+g(T.h)  (la) 

S° 

where  g  =  G/N  and  the  scale  factor  s  e  (N/N1)*^.  In  Eq  (la)  h'  and  g 

are  not  yet  fully  determined.  Now  take  the  thermodynamic  limit  N,N'  -*». 

We  now  make  our  main  assumption,  namely  that  near  the  transition  (e.g. 

for  small  h)  it  Is  possible  to  pick  h'  so  that  h'  and  g  are  (sufficiently) 

7  8 

regular  functions  of  h  (when  s«»)  *  . 

Now  for  values  of  T  for  which  a  first-order  transition  exists,  m 

q 

will  be  a  discontinuous  function  of  h  at  some  point  which  we  may  take 
as  h  =  0.  Differentiating  Eq  (1)  with  respect  to  h  at  h  *  0  +  e  we  find 

m  =  (lF‘)  (2) 

where  Am  Is  the  discontinuity  of  the  magnetization  at  temperature  T. 

No  term  In  dg/3h  appears  In  Eq.  (2)  since  g  Is  assumed  regular  at  h  =  0. 

Hence,  at  a  first  order  transition  in  an  infinite  system. 


3h' 

TFT 


(3) 


This  Is  essentially*®  the  result,  of  Nlenhuis  and  Nauenberg^ .  When  f 
Is  singular  but  m  goes  to  zero  as  some  power  hp  (p>0,  pf  Integer),  as  It 
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does  at  a  higher  order  phase  transition,  an  extension  of  this  argument 
shows  that  the  power  d  In  Eq(3)  Is  replaced  by  y  =  d/p+1  <  d.  Thus  the 
magnetic  field  Is  maximally  renormalized  at  a  first  order  phase  transi¬ 
tion^. 

It  has  .been  sbown  iby  Nightingale? -and  others 1?  that  very  accurate 
results  for  the  transition  temperature  and  critical  exponents  at  a  higher 
order  phase  transition  may  be  obtained  from  finite  system  results,  by 
properly  comparing  two  systems  of  nearby  sizes!*.  The  argument  that 
follows  extends  Nightingale's  method  to  first  order  transitions. 

Consider  Eq  (1)  for  finite  size  systems.  Differentiating  with 
respect  to  the  field,  we  have 


(h)  =  -jj  |j~  mN,  (h ' )  +  mN  (h).  (4) 

S 

The  tilde  in  Eq  (4)  and  below  refers  to  derivatives  of  g.  We  have 

ah' 

suppressed  T  and  also  the  dependence  of  on  N  and  N'.  In  general, 
for  finite  N  and  N',  all  terms  in  Eq  (4)  are  analytic.  To  obtain  a 
useful  result,  we  differentiate  with  respect  to  h  once  again.  This  gives 

XN  (h)  =  ^  <!f)2  xN'  (h'}  +id^r"V  (h,)  +  XN  (h)’  (5) 

2  2 

where  the  magnetic  susceptibility  x  =  3m/3h  *  -3  f/3h  .  Now  mN,  is 

(for  most  systems)  bounded  in  the  thermodynamic  limit  and  by  assumption 
-  2  2 

x  and  3  h'/3h  remain  finite.  At  a  first  order  phase  transition, 

as  N  and  N'  -*«,x  will  diverge  and  taking  leading  terms  we  have 


_  1  /3h 1 \ 2  pd 

XN  ~3  (3TP  xN'  -  S  XN' 
XN  *  F'  XN' 


(6) 
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Note  that  If  we  are  not  at  a  first  order  transition,  either  we  have  a 
higher  order  transition  for  which  dh ' /ah  =  S^,  y  <  d,  which  results  In 
a  lower  power18  of  N/N*  In  Eq  (6),  or  else  does  not  diverge  at  all. 

So  we  may  write 

XN  1  (Jr)  XH'.  for  N'  <  N,  i  (7) 

where  the  equal  sign  holds  for  first  order  transitions  only18,17*18.  Thus 
Eq  (7)  provides  a  criterion  for  locating  such  transitions  by  comparing 
finite  size  results.  It  also  provides  a  definition  of  a  first  order  transi¬ 
tion  as  a  thermodynamic  state  for  which  the  susceptibility  diverges  propor¬ 
tional  to  the  volume  of  the  system.  Note  that  it  implies  that  in  such  a  state, 

a  magnetic  field  of  0(1 /N)  is  sufficient  to  cause  a  finite  magnetization  in 

1 9 

the  thermodynamic  limit  . 

To  test  the  usefulness  of  our  criterion  (Eq.  7),  we  applied  it  to  the 
two  dimensional  ferromagnetic  Ising  model  on  a  lattice  with  square  unit  mesh 
with  nearest  neighbor  spin  coupling  J  *  1  and  periodic  boundary  conditions17. 
This  system  has  a  line  of  first  order  phase  transitions  for  (external  field) 
h  =  0  and  T  <  Tc  =  2.269  and  a  second  order  transition  at  T  =  Tc-  Using 
the  transfer  matrix  method,  we  determined  for  n  x  n  lattices  with  n 

p 

between  4  and  8  (note  N  =  n  ).  This  was  done  by  calculating  mN  for  sufficiently 
small  h,  so  that  the  susceptibility  could  be  obtained  from  xN  =  mN/h. 

For  T  values  well  above  Tc,  a  plot  of  x^  vs.  N  has  a  slope  that  decreases 
continually  as  N  increases.  For  T  well  below  Tc,  x^  becomes  linear  with 
increasing  N,  and  the  slope  increases  with  decreasing  T.  The  existence  of  a 
non-zero  slope  as  N-*»  indicates  (by  Eq  7)  that  one  is  in  a  first-order 
transition  region.  As  one  moves  out  of  such  a  region,  the  slope  will  drop 
more  or  less  abruptly  to  zero.  Finite  size  effects  smooth  out  this  drop. 

Hence  the  boundary  of  the  region  may  be  estimated  by  finding  the  Inflection 
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3Xu  3XN 

point  of  Letting  T*  be  the  temperature  for  which  (jr-)  ■  0  we 

find,  for  the  Islng  model,  T|  =  2.40,  T*  *  2.39,  Tf  *  2.37,  T*  *  2.36, 

20 

compared  to  T*  *  2.27  for  the  Infinite  system  . 

A  somewhat  more  accurate  result  may  be  obtained  by  calculating  xn  for 

an  n  x  »  strip  of  spins.  Here,  at  a  first  order  phase  transition,  xn  may 

grow  very  rapidly  with  n  since  the  system  Is  unbounded  In  one  direction.  In 

such  a  region  the  system  forms  many  domains  of  predominantly  up  or  predominantly 

down  spins.  For  a  given  value  of  n,  at  low  enough  temperatures,  the  domains 

will  form  a  one  dimensional  array  with  average  domain  length  (mean  ordered 

length)  iQ  -  exp  I&no),  where  o  is  the  free  energy  per  spin  of  a  domain  wall  . 

22 

Here  one  expects  xn  =  n  .  When  one  is  not  at  a  first  order  transition, 
xn  will  approach  a  constant  value  (or  diverge  less  rapidly)  as  n  -»®.  The  change¬ 
over  to  a  divergent  susceptibility  whose  rate  of  divergence  increases  rapidly 

with  temperature  should  allow  a  more  accurate  location  of  the  first  order 

2 

transition  than  for  the  n  x  n  lattice,  for  which  x  can  diverge  as  n  at  most. 

For  the  Ising  model  described  on  an  n  x  «  lattice,  we  find  that  at  low  tempera¬ 
tures  (xn/n)^n  apparently  tends  to  a  non-zero  constant  value  with  increasing 
n.  Plotting  Axn  =  xn/n  -  xn»/(n’)  »  n*  =  n-1,  one  finds  a  change  from 
negative  Ax  values  to  positive  Ax  values  at  the  point  where  xn  begins  to 
diverge  more  rapidly  than  n  .  Defining  7^  as  the  temperature  at  which  Ax  =  0, 
we  find  Tg  ■  2.209,  Tg  »  2.225,  Tj0  *  2.233,  compared  to  T'  *  2.269  for  the 
fully  Infinite  lattice.  In  many  systems  It  Is  not  possible  to  calculate 
with  one  dimension  Infinite.  However  one  can  probably  use  a  variant  of  this 
method  with  one  dimension  longer  than  the  other  (i.e.  a  rectangular  or  parallel - 
piped  shaped  lattice)  to  good  advantage. 

It  Is  clear  that  our  argument  may  be  extended  to  systems  other  than 
23  24 

Islng  models  *  .  In  general,  at  a  first-order  phase  transition  there  will 
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be  an  order  parameter  P  which  Is  a  discontinuous  function  of  a  conjugate 
thermodynamic  field  £,  where 

P  =  df/d£. 

The  discontinuity  will  occur  In  the  Infinite  system.  jOne  expects  that 

near  the  transition  f  has  the  shape  of  a  roof  top  wit^i  £  =  0  defining  the 

ridge  line.  Then  one  finds,  for  the  finite  system,  that  Eq  (6)  will  hold 

2  2 

with  x  replaced  by  3  f/3  £.  If  £  depends  on  the  temperature,  the  entropy 

26 

will  be  discontinuous  -  I.e.  there  will  be  a  latent  heat  of  transition  . 

In  this  case  one  expects  that  Eq  (3)  applies  with  h'  and  h  replaced  by  T'  and 

T,  respectively,  and  Eq  (6)  applies  for  the  specific  heat  C  ^  i.e.  C  diverges 

proportional  to  the  volume.  Eq  (3)  then  implies  that  the  specific  heat 

critical  exponent  a+1  at  a  first  order  transition  with  latent  heat.  This 

27 

conclusion  has  also  been  arrived  at  by  Imry  by  means  of  a  plausibility 
argument  for  the  finite  size  broadening  of  a  first  order  phase  transition. 

We  would  like  to  thank  G.  Akinci  for  calculating  the  results  for  the  n  x  n 
lattice  reported  above  and  for  some  very  insightful  remarks.  It  is  also  a 
pleasure  to  acknowledge  useful  discussions  with  R.  B.  Griffiths,  M.  E.  Fisher, 
J.  L.  Lebowitz,  E.  Lieb,  A.  N.  Berker  and  some  helpful  comments  from  one  of 
the  referees. 
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FOOTNOTES 

1.  B.  Nienhuis  and  M.  Nauenberg,  Phys.  Rev.  Letters  35,  477  (1975). 

2.  M.P.  Nightingale,  Physlca  8 3A,  561  (1976);  Phys.  Lett,  59A,  486 
(1977);  Proc.  Kon.  Ned.  Akad.  Wet.  82  series  B,  235  (1979). 

3.  R.H.  Swend sen,  Phys.  Rev.  B20,  2080  (1979), 

2 

4.  The  result  is  al9o  true  as  an  upper  bound  (i.e.  »  see  below 

for  notation)  for  lattice  systems  in  which  the  order  parameter  is 
the  expectation  value  of  an  operator  that  commutes  with  the  rest  of 
the  Hamiltonian  (R.B.  Griffiths,  Phys.  Rev.  1 52,  240  (1966)). 

5.  T.D.  Schultz,  D.C.  Mattis  and  E.  Lieb,  Rev.  Mod.  Phys.  38,  856  (1964) 
point  out  that  the  calculation  of  the  magnetization  of  the  two- 
dimensional  Ising  model  by  C.N.  Yang,  Phys.  Rev.  85,  808  (1952) 

(and  also  that  of  E.W.  Montroll,  R.B.  Potts,  and  J.C.  Ward,  J.  Math. 
Phys.  4,  308  (1963))  is  in  fact  a  calculation  of  /x/N  for  N  -*». 

Thus  x  *  N  when  T  <  Tc-  This  holds  for  periodic  boundary  conditions, 
for  which  the  zerofield  susceptibility  per  spin  xN  (see  Eq  (5)  and 
below)  may,  since  the  magnetic  field  h  is  exactly  zero,  be  written 
as  a  sum  of  N  two-spin  correlation  functions  <SiS.>,  since  the 
magnetization  vanishes  identically.  This  result  may  be  extended 
to  any  Ising  model  In  two  dimensions  with  purely  even  ferromagnetic 
interactions  including  nearest  neighbor  exchange  J  and  periodic 
boundary  conditions  by  using  the  GKS  Inequalities  (R.B.  Griffiths, 
in  Phase  Transitions  and  Critical  Phenomena,  C.  Domb  and  M.S.  Green, 
eds..  Academic  Press  (1972),  Vol.  1,  pp.  7-109)  since  each  <S^Sj> 

Is  bounded  below  by  <S^Sj>  with  J  only.  In  three  dimensions, 
one  can  embed  a  suitably  deformed  (to  preserve  the  boundary 
conditions)  plane  including  the  points  1  and  j  for  each  term  in  the 
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sum.  Then  for  T  <  T  one  has  x  ^N(m*)  where  *  refers  to  the  2d 

w 

Ising  model  with  interaction  J  only.  Higher  spin  models  may  be 
treated  using  the  transformation  of  Griffiths  (op.  cit.  and  J.  Math. 
Phys.  1£,  1559  (1969)).  Note  that  the  result  x  “  N  does  not  follow 
from  <S^Sj>  -♦in*  as  N  ♦  -  since  the  range  of  distances  jli-jl  over  - 
which  this  occurs  may  not  grow  sufficiently  rapidly  with  N.  However, 
the  power  N  is  clearly  maximal  whenever  x  is  a  sum  of  (bounded) 
correlation  functions  (cf.  footnote  4  and  12)  since  there  are  only 
N  terms  in  the  sum. 

6.  See,  for  instance,  Th.  Niemeijer  and  J.M.J.  Van  Leeuwen  in  Phase 
Transitions  and  Critical  Phenomena  Vol.  6,  C.  Domb  and  M.S.  Green, 
eds..  Academic,  New  York  (1976). 

7.  For  most  applications  of  our  criterion  (see  below)  it  is  probably 
sufficient  to  let  N  -  N'  be  fixed  or  otherwise  such  that  as  N  -*•  « 
our  assumptions  need  only  hold  for  s  ♦  1. 

2  2  2  2 

8.  "Sufficiently  regular"  means  that  we  require  3  h'/3h  and  3  g/3h 
to  be  finite  at  the  transition.  Thus  singularities  appearing  in 
higher  derivatives  of  h*  or  g  or  essential  singularities  are  com¬ 
patible  with  our  assumption. 

9.  This  condition  defines  a  first  order  phase  transition  for  our  purposes. 
Note  that  if,  for  small  h,  f  has  the  simple  form  f(T,h)  =  A(T)  + 
(am(T)/2)  |h|  +  B(T)h,  the  conditions  on  h'  and  g  are  satisfied  by 

h'  *  s^h,  g  ■  A(T)  (l-s-<*).  This  form  for  f  is  of  course  just  one 
possibility  and  by  no  means  unique. 

10.  Since  our  argument  rests  on  a  scaling  assumption  (E^  (la)  and  ff) 
it  is  not  necessary  to  Invoke  a  complete  microscopic  renormalization 
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group  mechanism  as  in  ref.  1.  Note  also  that  our  assumption  that 
only  one  parameter  (in  this  case  h)  is  renormalized  follows  that  of 
the  phenomenological  renormalization  group  treatment  of  Nightingale 
j(ref.  2)  for  higher  order  transitions. 

11.  jNote  that  Eqs.  (la)  and  (3)  are  consistent  with  assuming  that  the 
singular  part  of  f  is  a  homogeneous  function  of  h  with  degree  1  at 
a  first  order  phase  transition.  This  assumption,  in  turn,  is  a 
special  case  of  Widom's  scaling  equation  of  state  given  in  J.  Chem. 
Phys.  43,  3898  (1965). 

12.  In  renormalization  group  theory,  y  <  d  is  a  necessary  condition  for 
a  physical  spin  coupling  (one  that  does  not  grow  with  separation). 

In  scaling  theory,  y  >  d  would  imply  that  a  first  derivative  of  f 

is  infinite  at  the  transition,  which  is  unphysical.  (See  also  foot¬ 
note  5). 

13.  W.  Kinzel  and  M.  Schick,  Phys.  Rev.  B,  to  appear;  L.  Sneddon,  J. 
Phys.  C.  12,  3051  (1972);  Z.  Racz,  Phys.  Rev.  B21,  4012  (1980). 

14.  As  explained  in  ref.  2,  one  may  regard  the  change  in  size  of  the 
system  (from  N  to  N1  spins)  as  a  renormalization  group  transformation 
with  scale  change  s.  Hence,  one  can  in  principle  derive  a  scaling 
form  for  f  at  a  first  order  transition  by  the  standard  argument 

(see  ref.  4). 

15.  At  a  second  order  transition,  Sd  is  replaced  by  S2-n  in  Eq  (6). 

See  M.E.  Fisher  in  Critical  Phenomena,  M.S.  Green,  Ed.,  Academic, 

NY  (1971). 

16.  If  f  has  the  simple  form  given  in  footnote  7,  its  singular  part 
Is  «  |h|.  It  is  tempting  to  suppose  that  for  finite  N  this  is 
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simply  replaced  by  /h  +  hfi^  (N)  .  Eq  (6)  then  implies  hg^  ® 

1/N,  or  more  completely  heff=4m/(2x0N)»  where  x  =  XQ  (T)  N. 

For  a  transition  with  latent  heat,  this  functional  form  leads  to 
a  finite  size  rounding  aTc  of  the  transition  temperature.  Making 
I  the  identification  Tg^  =  aTc,  we  find  aTc  «  1/N,  in  agreement 
!  with  the  result  of  Imry  (ref.  27).  His  argument  further  implies 
that  T  ^  ^  kTc/(ASN)  where  AS  is  the  entropy  discontinuity. 

17.  M.E.  Fisher  has  remarked  (private  communication)  that  one  must  be 
extremely  careful  in  choosing  the  proper  boundary  conditions  in 
order  that  x  diverge  *  N.  We  are  indebted  to  one  of  the  referees 
for  pointing  out  that  the  intuitive  requirement  for  this  is  that 
both  phases  be  present  with  comparable  probabilities  -  otherwise 
the  variance  of  the  order  parameter  will  be  greatly  reduced.  If 

for  instance,  one  has  an  Ising  model  at  h=0  with  all  boundary  spins  up, 
in  general  m^  >  0  and  x  will  not  diverge.  The  work  of  J.L.  Lebowitz, 

J.  Stat.  Rhys.  16,  463  (1967)  is  interesting  in  this  regard. 

18.  The  result  that  x  diverges  «  N  as  well  as  several  other  general 
features  of  first  order  phase  transitions  have  recently  been  obtained 
by  A.N.  Berker  and  M.E.  Fisher  (unpublished). 

19.  This  point  has  been  demonstrated  for  the  Kac  model  (a  one  dimensional 
Ising  model  with  long  range  interactions)  and  the  analagous  point  for 
the  free  Bose  gas  by  Y.  Imry,  Ann.  Phys.  51^,  1  (1969). 

20.  Finding  the  inflection  point  of  3x^/9N  with  respect  to  6  =  1/kgT  instead 

* 

of  T  increases  the  error  in  Tg  by  50%. 

21.  L.  Onsager,  Phys.  Rev.  65,  117  (1944). 

22.  It  may  seem  somewhat  Inconsistent  to  use  a  semi-infinite  system,  for  which 
our  scaling  argument  does  not  apply,  in  order  to  get  (Improved)  numerical 
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results.  However,  the  same  procedure  is  employed  in  the  "phenomeno¬ 
logical  renormalization  group"  theory  of  Nightingale  (Refs.  2  and  13). 

In  this  case  one  calculates  the  correlation  length  £  as  1 /I n ( /Xq) . 
where  Xg  is  the  largest  and  X-j  the  second  largest  eigenvalue  of  the 
transfer  matrix.  This  formula  gives  the  correlation  length  correctly 
for  an  n  x  «  strip  in  a  disordered  region,  and  diverges  «  n  a(  a  second 
order  phase  transition,  in  accord  with  scaling  theory.  However  at  a 
first  order  transition  (e.g.  for  h  ■  0,  T  <  Tc  in  the  Ising  model), 

€  calculated  according  to  this  formula  diverges  exponentially  with  n. 

This  is  because  it  is  then  a  measure  of  the  size  of  the  ordered  domains 
referred  to  above,  closely  related  to  the  mean  ordered  length  Hq.  Were 
one  to  calculate  the  correlation  length  on  a  finite  lattice  (e.g.  an  n  x  n 
square)  it  would  in  general  diverge  less  rapidly  (or  not  at  all)  in  the 
first  order  region,  and  we  expect  this  would  reduce  the  accuracy  with 
which  the  second  order  transition  can  be  located,  similarly  to  what 
we  find  herein  for  the  first  order  case. 

23.  L.K.  Runnels,  J.R.  Craig  and  H.R.  Strieffer,  J.  Chem.  Phys.  54,  2004  (1971) 
conclude  that  there  is  a  first-order  phase  transition  in  a  certain 
triangular  lattice  gas  with  hard  cores  on  an  n  x  ~  lattice.  Their 
results  indicate  that  at  the  transition  the  compressibility  diverges 
faster  than  n2  (as  n2’^)  for  a  certain  finite  range  of  values  of  n, 
consistent  with  our  results. 

24.  However,  M.E.  Fisher  (private  communication)  suggests  that  there  may 
be  certain  systems  (e.g.  with  vector  spins)  for  which  Am  f  0  but  x 
does  not  diverge  «  N  (e.g.  when  the  long-long  range  order  vanishes). 

Thus  although  the  derivation  of  Eq  (6)  Is  very  general,  its  ultimate 
range  of  validity  remains  to  be  established. 
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25.  Note  that  the  definition  of  P  may  not  be  obvious  In  all  cases  e.g.  if 
there  Is  no  obvious  symmetry.  This  might  pose  a  problem  In  the  correct 
choice  of  boundary  conditions  (see  Ref.  17). 

26.  A,simple  spin  model  illustrating  this  point  is  provided  by  R.B. 
Griffiths,  Physica  33^  689  (1967). 

27.  yJ  Irary,  Phys.  Rev.  B21^,  2042  (1980).  It  is  easy  to  show  that  this 
work  actually  implies  c  =  (as)  N/k,  in  agreement  with  Eq  (6),  where 
as  Is  the  entropy  discontinuity. 


